%% ——— 模型微分方程 ———
function dS = f(S, alpha, beta, xi, mu, mu_s)
    x = S(1); y = S(2); z = S(3);
    dx = -alpha*( y - x - 0.5*(y-x)^2 - (y-x)^3/3 + beta*sin(z) );
    dy = mu_s - xi*y + ( y - x - 0.5*(y-x)^2 - (y-x)^3/3 );
    dz = mu * x;
    dS = [dx; dy; dz];
end